Polarizability of molecular chains: does one need exact exchange? 
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Standard density functional approximations greatly over-estimate the static polarizability of long- 
chain polymers, but Hartree-Fock or exact exchange calculations do not. Simple self-interaction 
corrected (SIC) approximations can be even better than exact exchange, while their computational 
cost can scale only linearly with the number of occupied orbitals. 
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Ground-state Kohn-Sham (KS) density functional the- 
ory (DFT) has become extraordinarily popular for solv- 
ing electronic structure problems in solid-state physics, 
quantum chemistry and materials science [l[. The ac- 
curacy of modern generalized gradient approximations 
(GGAs) and hybrid functionals has proven sufficient for 
many applications, often with surprisingly small errors. 
Bond dissociation energies, geometries, phonons, etc., are 
now routinely calculated with errors of 10-20%. 

But local and gradient-corrected functionals overesti- 
mate massively the static polarizability and hyperpolar- 
izability of molecular chains, especially conjugated poly- 
mers. This failure has been the subject of many studies 
over the last decade @, i, S U H 0, E & studies 
which highlight the important role played by the response 
field originating from the exchange-correlation (XC) po- 
tential. The exact induced XC field counteracts the ap- 
plied external field, keeping the polarization low. In 
the local (or gradient-corrected) density approximation 
(LDA), this field erroneously points in the same direc- 
tion as the applied field [!, 0, H[ . Such failures of stan- 
dard functionals appear in other contexts, such as trans- 
port through single molecules llj, or the polarizability 
of large molecules. 

In contrast, these effects are easily captured within 
standard wavefunction theory. In particular, Hartree- 
Fock (HF) theory does not greatly overestimate the po- 
larizabilities and provides a good starting point for more 
accurate wavefunction treatments, such as Moller-Plesset 
(MP) perturbation theory. Thus exact exchange (EXX) 
DFT, the KS-DFT method for minimizing the HF energy 
while retaining a single multiplicative potential, provides 
a promising alternative and indeed has been found to give 
results very similar to HF 0, H, E3] • This improvement 
can be attributed to the orbital-dependence of EXX, and 
the lack of self-interaction error ||, i.e. EXX is exact for 
one electron, unlike LDA or GGA. 

However, EXX is only one among many possible self- 
interaction free functionals that one may construct. In 
fact any GGA can be corrected to become self-interaction 
free (self-interaction corrected - SIC) by direct subtrac- 
tion of the XC functional evaluated on each of the indi- 
vidual orbitals [l2|. While this can be performed for ei- 
ther LDA or GGA, only LDA has significantly improved 
energetics from this procedure, but many investigators 



are searching for useful methods to correct GGA's for 
self- interaction jl3j |. More importantly, EXX includes a 
sum over all unoccupied orbitals, while SIC functionals 
use only occupied ones, i.e., EXX can often be substan- 
tially more expensive computationally. So the question 
then becomes: does one really need EXX, or will any 
self-interaction free functional perform equally well ? 

We perform SIC calculations for the polarizabilities of 
hydrogenic chains using LDA and GGA. Our SIC po- 
tential is constructed using the optimized effective po- 
tential (OEP) framework within the Krieger-Li-Iafrate 
(KLI) approximation [ijj]. Using results from accurate 
wave-function methods as a benchmark, we find that 
the polarizabilities calculated with KLI-SIC are in bet- 
ter agreement than those obtained with KLI exact ex- 
change (X-KLI), with the remaining error attributed to 
the KLI approximation [Tol . [To] ]. This is an important 
result since KLI-SIC scales only linearly with the num- 
ber of occupied Kohn-Sham (KS) orbitals, as compared 
with the quadratic scaling of X-KLI. Thus SIC becomes a 
valuable scheme for evaluating the static polarizability of 
polymers with large unit cells, and in other applications 
where these effects may be important [llj]. 

We start with a brief description of the SIC method 
used in this work. In DFT [la] , the total energy func- 
tional E[p\p^] (p a is the spin a =T,J. density, p = 
J2<j P a ) can be written as 

E[p\pt] =T s [p] + J d 3 rp(r)v(r) + U{p} + E xc [p\pt}, 

(1) 

with Ts the kinetic energy of the non-interacting KS or- 
bitals, v(r) the external potential, U the Hartree energy, 
and E xc the XC energy. For any GGA 
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where p° — |"0n| 2 is the density of the n-th KS orbital. 
Levy's minimization [I?) leads to a set of single particle 
KS-likc equations for with corresponding eigenvalues 
e^ ,SIC and occupation numbers "p a n (p a = J2 n PnPn) 
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The effective potential v° s (r) is now KS-orbital de- 
pendent and cannot be classified as a standard multi- 
plicative KS-potential. For instance it is ambiguously 
denned for unoccupied KS-orbitals since the SIC is only 
defined for the occupied ones. The solution of equation 
([3|) has followed several approaches. The simplest one 
is to solve it directly under a normalization constraint 
with the resulting non-orthogonal orbitals undergoing an 
orthogonalization procedure [l2|. This scheme however 
is not free of complications, since E^P is not invariant 
under a unitary transformations of the occupied {ip^}- 
The solution [181 ] is then to work with an auxiliary set 
of localized orbitals used for constructing v^ s n (r) 

and related to {V'nl by a unitary transformation chosen 
so as to minimize £ X S C IC 

A convenient alternative is offered by the OEP method 
[l9l ] where the orbital-dependent v° s n (r) is recast into 
a local multiplicative orbital-independent potential. In 
this way the SIC problem can then be solved as a nor- 
mal KS problem. However the construction of an OEP 
is computationally demanding. Here we adopt the KLI 
approximation [l4| , which is practically easy to construct 
and retains most of the advantages of the full OEP, be- 
ing therefore well suited to the SIC problem [20]. The 
orbital-independent KLI effective SIC Kohn-Sham poten- 
tial takes the form 

t£ s (r) = v(r) + vn(r) + <' GGA (r) + <- SIC (r), (4) 

where «h and v^ GGA are the Hartree and GGA-XC po- 
tentials. We define the SIC potentials by 

<' SIC (r) = -^' x G c GA [/5« ; 0](r), (5) 

where t>Hxc is the sum of the Hartree and GGA potentials, 
and their Slater average as 

N° 

<c SIC (r) = E/«( r K' SIC ( r )' ( 6 ) 

71=1 

where /£(r) = p^(r)/p CT (r) is the auxiliary orbital den- 
sity as a fraction of the spin-density of its spin. Then 

N" 

<> SIC (r) = <f IC (r) + Y, £(r)[At# SIC - C*]. (7) 

n=l 

The orbital densities p^ in Eqs. ©-([7]) are calculated 
from the auxiliary set of localized orbitals {4>f} instead 
of the canonical {tpf }. Both sets of orbitals give the same 
total density p" . The orbital shift terms Av^ being a 
constant), are obtained by solving 

N " —cr 

^ ] (5nrn — fnm)^ v n ~ w xc,m ~ u m J (8) 

n=l 

where \i m = / dr j 5^ l (r)/i <7 (r). Furthermore, 

N" _ N" 

£7L = i, E«^ c -<i SIC )=o, (9) 

71=1 71=1 



so that the linear system in Eq.© is of rank (N^-l) 
and must be solved using a least squares approach. The 
constant C a is set to Ad^'q^q (HOMO is the highest 
occupied molecular orbital). We have implemented this 
scheme in the DFT code SIESTA 

To facilitate easy comparison with pre vious quantum 
chemistry and EXX-DFT calculations [10| , we chose as a 
test system the widely studied linear hydrogen chains H„ 
made up of n H atoms with alternating H-H distances of 
2 and 3 a,Q. An optimized atomic orbital basis set consist- 
ing of double zeta polarized s and triple zeta polarized 
p functions is employed. The Pipek-Mezey localization 
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FIG. 1: The Pipek-Mezey localized molecular orbitals {4>i} 
(left) and canonical Kohn-Sham orbitals {ipl} (right) for the 
occupied states of the Hg molecule. 

scheme [22| which minimizes the number of atoms over 
which a given molecular orbital is delocalized, is used to 
transform the canonical KS-orbitals {ip^} into the local- 
ized set {4>Z}- As an example both the {c^} and the 
{tpn} °f ^6 are presented in figure [TJ Finally the static 
polarizability a — 5p z /8F z is calculated numerically us- 
ing finite differences. 

Table U contains the central results of this work. The 
first column contains highly accurate (MP4) quantum 
chemical results, which we take as exact. The next two 
columns show LDA and PBE results [HI], demonstrat- 
ing the LDA overestimate (by about 100% for H12) of a, 
an overestimate that is only slightly reduced by GGA. 
We checked several GGA's, and they all had the same 
features. Small finite systems, such as atoms, are well- 
known to overpolarize in LDA and GGA, because their 
XC potentials are too shallow, and do not decay with 
the correct asymptotic form, — 1/r. But this is a distinct 
effect from that discussed here, which is due to the re- 
sponse to the electric field inside the molecule. Ours is a 
bulk effect, not depending on the end points. We checked 
this and found that the LB94 functional [24[ , specifically 
designed to reproduce the correct asymptotic behavior, 
does not yield any better results than the other GGA's 
and in fact worsens the LDA a. 

The next 3 columns list results of different types of 
calculations using the Fock integral, and no correlation. 
We note that HF slightly overestimates the polarizabil- 
ity, but by less than 10%. An exact OEP treatment of 
the same functional (X-OEP) yields essentially the same 
numbers. Our KLI scheme, applied to the same func- 
tional (X-KLI), makes a noticeable overestimate, but the 
error remains less than 20%, compared to OEP. 

Now we focus on the SIC results. Regardless of the 
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TABLE I: Calculated polarizability a of H n chains obtained by using the KLI-SIC method with different XC functionals. The 
subscript X indicates that correlation has been dropped from the XC potential. These are compared with MP4, HF and exact 
exchange DFT (EXX) results from reference [10|]. 





exact 


GGA 




EXX 






KLI-SIC 






X-only 




Hjv 


MP4 


LDA PBE 


HF 


X-OEP 


X-KLI 


LDA 


PBE 


LDA-OEP a fl5] 


LDA X 


PBE X 


SIC-LDAx 


SIC-PBE X 


H 4 


29.5 


37.26 35.62 


32.0 


32.2 


33.11 


33.38 


33.14 


32.5 


38.90 


36.51 


33.37 


33.10 


H 6 


51.6 


73.10 69.35 


56.4 


56.6 


60.64 


58.56 


58.07 


54.6 


76.16 


70.45 


58.84 


57.63 


H 8 


75.9 


116.58 109.74 


82.3 


84.2 


91.56 


86.94 


86.48 


76.4 


121.64 


110.84 


87.08 


84.53 


Hio 




166.36 155.31 






124.87 


117.28 


116.16 




173.89 


156.45 


116.77 


113.14 


H12 


126.9 


220.55 204.53 


137.6 


138.1 


159.27 


147.96 


145.98 


126.8 


231.25 


205.51 


147.19 


141.90 



"Estimated by subtracting the difference between X-KLI and X- 
OEP. 



GGA functional used, the KLI-SIC polarizabilities show 
a drastic improvement with respect to those obtained 
with pure GGA. Furthermore, PBE results are very close 
to LDA results. There remains about a 20% overestimate 
for H12, but this is probably due to the KLI approxi- 
mation. Correcting the KLI-SIC-LDA values, using the 
difference between KLI and OEP for EXX as an error 
estimate, suggests that OEP-SIC rcsults(C/. Korzdorfer 
et al.^\j^) might be in good agreement with MP4. These 
corrected values are listed under LDA-OEP in Table [H 

We also checked if the inclusion of correlation was im- 
portant for the excellent KLI-SIC results, by running the 
calculations with correlation removed. For the pure func- 
tionals, LDA correlation slightly reduces the huge over- 
estimate, whereas PBE correlation has almost no effect. 
On the other hand, for the KLI-SIC results, LDA correla- 
tion has almost no effect, while PBE correlation corrects 
the PBEx result in the wrong direction. This is consis- 
tent with the general result that SIC-GGA includes some 
incorrect overcounting of gradient effects |l3| . 
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FIG. 2: AFxc(z) = Vi c (z) - Vxc for H 8 plotted against 
position (z) for (a) LDA, (b) SIC-LDA, (c) X-KLI, (e) X- 
KLI S and (f) SIC-LDA S . The applied external field is shown 
in (d) and the black dots indicate the position of the H atoms. 

The improved response obtained with orbital de- 



TABLE II: Static polarizability a of H n obtained from SIC 
and X-KLI where only the Slater term in the KLI potential is 
used (super-script S). Both XC and X-only results are shown 
for the SIC-LDA S and SIC-PBE S . 



Hjv 


SIC-LDA S 


SIC-PBE S 


SIC-LDA| 


sic-pbe! 


X-KLI S 


H 4 


35.37 


35.12 


36.25 


35.16 


35.78 


H 6 


67.74 


67.13 


68.92 


66.34 


69.17 


H 8 


105.91 


104.89 


107.57 


102.97 


108.72 


Hio 


148.64 


146.18 


150.86 


143.10 


152.90 


H12 


193.94 


190.47 


197.05 


185.25 


199.91 



pendent functionals is due the opposing XC field, al- 
ready demonstrated in the literature for OEP exact 
exchange Here we verify that the same happens with 
the KLI-SIC scheme. In figure[2]we plot AVxc defined as 
the difference between the XC potential with and with- 
out an applied electric field AVxc( z ) = Vxc( z ) ~ ^xc for 
LDA, SIC-LDA and X-KLI. We denote by a superscript S 
results obtained by including only the Slater average po- 
tential in Eq. ([7]), dropping the constant terms. Clearly 
both contributions are significant in the final XC poten- 
tial. The Slater-only polarizabilities are shown in table 
ILTl and reflect this. 

By comparing the values of a in the tables Q] and |TT] 
one may conclude that the Slater average already con- 
tains important corrections, but the bulk of the effect is 
contained in the orbital shift term. For example if one 
considers the a calculated with LDA for H12 the polariz- 
ability is 220.5 for LDA, 193.9 for SIC-LDA S and 147.9 
for SIC-LDA. It is also interesting to note that even at the 
level of the Slater average, SIC performs better than X- 
KLI. This can be understood by looking at the XC poten- 
tial for SIC-PBE S and X-KLI S (Fig. |3]) when no external 
field is applied. The SIC-PBE potential exhibits higher 
peaks in the inter-molecular space between H2 units than 
X-KLI. This explains the quantitative difference in a be- 
tween the two cases. The improved performance of X- 
OEP over X-KLI can be attributed [10| to similar bar- 
riers in the inter-molecular region which however arise 
from the response part of the X-OEP potential. 

Finally we ask whether similar results can be obtained 
with atomic-like corrections, which have the effect of 
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making the KS potential deeper at the atomic sites. We 
investigated both the LDA+U (2f| and the atomic LDA- 
SIC (ASIC) methods [H,^ in this regard. The LDA+U 
results are extremely poor, as it provides polarizabilitics 
larger than even those obtained with simple LDA. The 
ASIC results are far more promising, being half-way be- 
tween those of GGA and of SIC (a=33.95, 63.67, 98.39, 
137.42, 178.87 for H„ respectively with n= 4, 6, 8, 10, 
12). In fact, if the ASIC corrections relative to pure 
LDA were double those we found, they would reproduce 
the HF results very accurately. 

The above can be understood by the way the atomic 
corrections are introduced. The H atoms are half-filled 
in the absence of an external field. Switching on the field 
produces a tiny charge transfer, which however is ampli- 
fied by the LDA+U potential. As a result the already 
wrong LDA response field gets amplified. The same does 
not happen with ASIC, which improves the response over 
LDA by virtue of the higher inter-molecular barriers. 

In conclusion we have demonstrated that SIC func- 
tionals at the KLI level, in general, perform better than 
X-KLI at computing the static polarizability of hydro- 
genic chains. This is an interesting result in view of the 
considerably less computational overheads involved with 



the SIC method with respect to EXX. Our work there- 
fore opens the prospect of using SIC for evaluating the 
electrical response of complex polymeric materials. 

This work is funded by the Science Foundation of Ire- 
land (SFI02/IN1/I175) and by the US Department of 
Energy (DE-FG02-01ER45928). 
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FIG. 3: XC potential plotted against distance (z) for SIC- 
PBE S (solid line) and X-KLI S . The grey inset plots the dif- 
ference between the two potentials. 
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